Intraspecific Trait Variation

Intro and packages

bhla bhla bhla

# data manipulation and visualization
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# correlation plot
library(corrplot)

## corrplot 0.94 loaded

# Functional diversity 
library(funspace); library(TPD)

## Loading required package: ks

Calling data

Measurements in the museum

data_specimens <- read_csv("flmnh_data/Museum Data Updated.csv") |>
  mutate(species_updated = paste(Genus, Species, sep = " "),
         body.mass = as.numeric(Mass))

## Rows: 284 Columns: 27
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (15): Meta-archipelago, Archipelago_group, Site, Genus, Species, Subspec...
## dbl (12): Area, DistanceM, Specimen ID, Wing.Length, Secondary.Length, Kipps...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `body.mass = as.numeric(Mass)`.
## Caused by warning:
## ! NAs introduced by coercion

species <- unique(data_specimens$species_updated)

AVONET raw data

We use the raw data (AVONET_Raw_Data sheet in the Excel file) in AVONET (Tobias et al. Ecol Letters 2022), to estimate a mean and SD of trait per species. Unfortunately, raw data of AVONET does not include the variation of Body mass, they only included the mean value in their summary (AVONET2_eBird sheet in the Excel file). We will pull the other sources of functional traits and estimate mean and SD of this trait afterwards. We also include the taxonomy information in AVONET (eBird).

avonet <- read_csv("FunctionalTraits/AVONET_Raw_Data.csv") |> 
  mutate(species_updated = ifelse(eBird.species.group == "Chlorophonia musica flavifrons",
                                  "Chlorophonia flavifrons",
                           ifelse(eBird.species.group == "Chlorophonia musica musica",
                                  "Chlorophonia musica",
                           ifelse(eBird.species.group == "Chlorophonia musica sclateri",
                                  "Chlorophonia sclateri",
                           Species2_eBird)))) |>
  dplyr::select(species_updated, Source, Specimen.number,
                Beak.Length_Culmen, Beak.Width, Beak.Depth,
                Tarsus.Length, `Hand-wing.Index`, Tail.Length) |>
  rename(HWI = `Hand-wing.Index`) |>
  filter(species_updated %in% species)

## Rows: 90371 Columns: 26
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (11): Avibase.ID, Species1_BirdLife, Species2_eBird, eBird.species.group...
## dbl (13): Data.type, Age, Beak.Length_Culmen, Beak.Length_Nares, Beak.Width,...
## lgl  (2): Locality, Country
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

This raw data does not have Mass values. We can extracted from the summarized AVONET2_eBird dataset.

avonet2_eBird <- read_csv("FunctionalTraits/AVONET2_eBird.csv") |>
  dplyr::select(Species2,
                Mass, `Hand-Wing.Index`, Tarsus.Length, Tail.Length, 
                Beak.Length_Culmen, Beak.Width, Beak.Depth) |>
  rename(body.mass = Mass) |>
  rename(HWI = `Hand-Wing.Index`)

## Rows: 10661 Columns: 31
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (13): Species2, Family2, Order2, Avibase.ID2, Mass.Source, Mass.Refs.Oth...
## dbl (18): Total.individuals, Female, Male, Unknown, Complete.measures, Beak....
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

avonet_raw <- avonet |>
  left_join(avonet2_eBird |> dplyr::select(Species2, body.mass), 
            join_by("species_updated" == "Species2"))

Finally, we added the archipelago group (manually) to these 383 specimens

avonet_islands <- read_csv("flmnh_data/AVONET_b_m_specimens_to_review.csv") |> 
  left_join(avonet_raw)

## New names:
## Rows: 383 Columns: 29
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (13): Avibase.ID, Species1_BirdLife, Species2_eBird, eBird.species.group... dbl
## (14): ...1, Data.type, Age, Beak.Length_Culmen, Beak.Length_Nares, Beak.... lgl
## (2): Locality, Country
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Joining with `by = join_by(Source, Specimen.number, Beak.Length_Culmen,
## Beak.Width, Beak.Depth, Tarsus.Length, Tail.Length, species_updated)`
## • `` -> `...1`

Combine data

The NA data will be replaced by the mean reported in AVONET

names(data_specimens)

##  [1] "Meta-archipelago"   "Archipelago_group"  "Site"              
##  [4] "Area"               "DistanceM"          "Genus"             
##  [7] "Species"            "Subspecies"         "Specimen ID"       
## [10] "Date"               "Data"               "Country"           
## [13] "Locality"           "Elevation"          "Other Info"        
## [16] "Sex"                "Mass"               "Wing.Length"       
## [19] "Secondary.Length"   "Kipps_D"            "HWI"               
## [22] "Tarsus.Length"      "Tail.Length"        "Beak.Width"        
## [25] "Beak.Depth"         "Beak.Length_Culmen" "Measurer"          
## [28] "species_updated"    "body.mass"

names(avonet_islands)

##  [1] "...1"                "Avibase.ID"          "Species1_BirdLife"  
##  [4] "Species2_eBird"      "eBird.species.group" "Species3_BirdTree"  
##  [7] "Data.type"           "Source"              "Specimen.number"    
## [10] "Sex"                 "Age"                 "Locality"           
## [13] "Country_WRI"         "Country"             "Beak.Length_Culmen" 
## [16] "Beak.Length_Nares"   "Beak.Width"          "Beak.Depth"         
## [19] "Tarsus.Length"       "Wing.Length"         "Kipps.Distance"     
## [22] "Secondary1"          "Hand-wing.Index"     "Tail.Length"        
## [25] "Measurer"            "Protocol"            "Publication"        
## [28] "species_updated"     "Archipelago_group"   "HWI"                
## [31] "body.mass"

traits_specimens <- data_specimens |> 
  dplyr::select(Archipelago_group, species_updated,
                body.mass, HWI, Tarsus.Length, Tail.Length, 
                Beak.Length_Culmen, Beak.Width, Beak.Depth)

traits_avonet <- avonet_islands |>
  dplyr::select(Archipelago_group, species_updated,
                body.mass, HWI, Tarsus.Length, Tail.Length, 
                Beak.Length_Culmen, Beak.Width, Beak.Depth)

data_traits <- traits_specimens |>
  rbind(traits_avonet) |>
  group_by(species_updated) |>
  mutate(body.mass = ifelse(is.na(body.mass), avonet2_eBird$body.mass,
                            body.mass),
         HWI = ifelse(is.na(HWI), avonet2_eBird$HWI,
                            HWI),
         Tarsus.Length = ifelse(is.na(Tarsus.Length), avonet2_eBird$Tarsus.Length,
                            Tarsus.Length),
         Tail.Length = ifelse(is.na(Tail.Length), avonet2_eBird$Tail.Length,
                            Tail.Length),
         Beak.Length_Culmen = ifelse(is.na(Beak.Length_Culmen), avonet2_eBird$Beak.Length_Culmen,
                            Beak.Length_Culmen),
         Beak.Width = ifelse(is.na(Beak.Width), avonet2_eBird$Beak.Width,
                            Beak.Width),
         Beak.Depth = ifelse(is.na(Beak.Depth), avonet2_eBird$Beak.Depth,
                            Beak.Depth))

Correct some names and add Id for each data, grouping by species

colnames(data_traits) <- c("Archipelago", "Species",
                           "mass", "h.w.i", "tars.l", "tail.l",
                           "beak.l", "beak.w", "beak.d")

data_traits <- data_traits |>
  group_by(Species, Archipelago) |>
  mutate(Sequence = row_number()) |>
  ungroup()

# a function for extract strings
id_specimen <- function(archipelago, species) {
    archipelago_part <- str_extract_all(archipelago, 
                                        "\\b\\w{2}") |>
      unlist() |>
      paste0(collapse = "")
    
    species_part <- str_extract_all(species, 
                                    "\\b\\w{3}") |>
      unlist() |>
      paste0(collapse = "")
    
    paste0(species_part, "_", archipelago_part)
}

data_traits <- data_traits |>
  mutate(Group = mapply(id_specimen, Archipelago, Species),
         ID = paste0(Group, "_", Sequence))

Functional Trait Space - PCA based

Chlorophonia-Euphonia

eup_traits <- data_traits |>
  filter(Species != "Coereba flaveola")

eup_traits_scl <- as.data.frame(scale(eup_traits[c(3:9)]))

corrplot(round(cor(eup_traits_scl),2), type="upper", order="hclust", 
         tl.col="black", tl.srt=45)

rownames(eup_traits_scl) <- eup_traits$ID

# Calculating the dimensionality of a functional space based on PCA
funspaceDim(eup_traits_scl)

## 
## Using eigendecomposition of correlation matrix.

## 
## Retain the first 3 components

## [1] 3

Run PCA of the scaled traits and building the functional trait spaces

pca.traits.eup <- princomp(eup_traits_scl, cor = T)

plot(pca.traits.eup)

And we can do a functional trait space analysis grouping by species and archipelago

# trait space by groups
trait.space.pca.eu.global <- funspace(x = pca.traits.eup,
                               n_divisions = 200)
# trait space by Archipelago
trait.space.pca.eu.arch <- funspace(x = pca.traits.eup,
                               group.vec = eup_traits$Archipelago,
                               n_divisions = 200)

# trait space by group of species by archipelago
trait.space.pca.eu.group<- funspace(x = pca.traits.eup,
                               group.vec = eup_traits$Group,
                               n_divisions = 200)

We can explore the results

summary(trait.space.pca.eu.global)

## 
## Functional space based on a PCA with 7 dimensions
## Dimensions 1 and 2 are considered in analyses
## 
## Loadings:
##        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mass    0.300  0.522  0.226  0.745  0.168  0.046  0.001
## h.w.i  -0.636 -0.229 -0.284  0.188  0.115  0.643 -0.018
## tars.l -0.360  0.566 -0.488  0.041 -0.557 -0.018  0.014
## tail.l -0.377  0.657 -0.168 -0.386  0.499 -0.023 -0.015
## beak.l -0.328 -0.335 -0.661  0.335  0.225 -0.424 -0.020
## beak.w -0.859 -0.038  0.387  0.078 -0.013 -0.136  0.294
## beak.d -0.825 -0.007  0.451  0.070 -0.089 -0.151 -0.284
## 
## Percentage of variance explained for each trait:
##        Comp.1 Comp.2 Overall_explained
## mass     9.03  27.27             36.30
## h.w.i   40.47   5.24             45.71
## tars.l  12.93  31.98             44.92
## tail.l  14.18  43.15             57.33
## beak.l  10.77  11.20             21.97
## beak.w  73.72   0.14             73.86
## beak.d  68.08   0.00             68.08
## 
## ------------------------------------------------------------
## 
## Functional diversity indicators:
## 
##  ---> For the global set of species:
## 
## Functional richness (99.9% probability threshold) = 47.26
## 
## Functional divergence  = 0.87

plot(x = trait.space.pca.eu.global, type = "global", 
     quant.plot = T, arrows = T, arrows.length = 0.75)

summary(trait.space.pca.eu.arch)

## 
## Functional space based on a PCA with 7 dimensions
## Dimensions 1 and 2 are considered in analyses
## 
## Loadings:
##        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mass    0.300  0.522  0.226  0.745  0.168  0.046  0.001
## h.w.i  -0.636 -0.229 -0.284  0.188  0.115  0.643 -0.018
## tars.l -0.360  0.566 -0.488  0.041 -0.557 -0.018  0.014
## tail.l -0.377  0.657 -0.168 -0.386  0.499 -0.023 -0.015
## beak.l -0.328 -0.335 -0.661  0.335  0.225 -0.424 -0.020
## beak.w -0.859 -0.038  0.387  0.078 -0.013 -0.136  0.294
## beak.d -0.825 -0.007  0.451  0.070 -0.089 -0.151 -0.284
## 
## Percentage of variance explained for each trait:
##        Comp.1 Comp.2 Overall_explained
## mass     9.03  27.27             36.30
## h.w.i   40.47   5.24             45.71
## tars.l  12.93  31.98             44.92
## tail.l  14.18  43.15             57.33
## beak.l  10.77  11.20             21.97
## beak.w  73.72   0.14             73.86
## beak.d  68.08   0.00             68.08
## 
## ------------------------------------------------------------
## 
## Functional diversity indicators:
## 
##  ---> For the global set of species:
## 
## Functional richness (99.9% probability threshold) = 47.26
## 
## Functional divergence  = 0.87
## 
##  ---> For each group:
## 
## Functional richness (99.9% probability threshold; Continental island) = 5.69 
## Functional richness (99.9% probability threshold; Greater Antilles) = 10.35 
## Functional richness (99.9% probability threshold; Lesser Antilles - Kalinago) = 3.48 
## Functional richness (99.9% probability threshold; Mainland) = 44.03
## 
## Functional divergence (Continental island) = 0.76 
## Functional divergence (Greater Antilles) = 0.83 
## Functional divergence (Lesser Antilles - Kalinago) = 0.39 
## Functional divergence (Mainland) = 0.84

plot(x = trait.space.pca.eu.arch, type = "groups", 
     quant.plot = T,globalContour = T, pnt = T, pnt.cex = 0.1)

summary(trait.space.pca.eu.group)

## 
## Functional space based on a PCA with 7 dimensions
## Dimensions 1 and 2 are considered in analyses
## 
## Loadings:
##        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mass    0.300  0.522  0.226  0.745  0.168  0.046  0.001
## h.w.i  -0.636 -0.229 -0.284  0.188  0.115  0.643 -0.018
## tars.l -0.360  0.566 -0.488  0.041 -0.557 -0.018  0.014
## tail.l -0.377  0.657 -0.168 -0.386  0.499 -0.023 -0.015
## beak.l -0.328 -0.335 -0.661  0.335  0.225 -0.424 -0.020
## beak.w -0.859 -0.038  0.387  0.078 -0.013 -0.136  0.294
## beak.d -0.825 -0.007  0.451  0.070 -0.089 -0.151 -0.284
## 
## Percentage of variance explained for each trait:
##        Comp.1 Comp.2 Overall_explained
## mass     9.03  27.27             36.30
## h.w.i   40.47   5.24             45.71
## tars.l  12.93  31.98             44.92
## tail.l  14.18  43.15             57.33
## beak.l  10.77  11.20             21.97
## beak.w  73.72   0.14             73.86
## beak.d  68.08   0.00             68.08
## 
## ------------------------------------------------------------
## 
## Functional diversity indicators:
## 
##  ---> For the global set of species:
## 
## Functional richness (99.9% probability threshold) = 47.26
## 
## Functional divergence  = 0.87
## 
##  ---> For each group:
## 
## Functional richness (99.9% probability threshold; Chlcya_Ma) = 9.65 
## Functional richness (99.9% probability threshold; Chlele_Ma) = 18.18 
## Functional richness (99.9% probability threshold; Chlfla_LeAnKa) = 3.48 
## Functional richness (99.9% probability threshold; Chlmus_GrAn) = 6.11 
## Functional richness (99.9% probability threshold; Chlocc_Ma) = 3.74 
## Functional richness (99.9% probability threshold; Chlscl_GrAn) = 8.5 
## Functional richness (99.9% probability threshold; Eupaff_Ma) = 7.93 
## Functional richness (99.9% probability threshold; Eupchl_Ma) = 7.95 
## Functional richness (99.9% probability threshold; Eupful_Ma) = 3.13 
## Functional richness (99.9% probability threshold; Eupgou_Ma) = 11.34 
## Functional richness (99.9% probability threshold; Euphir_Ma) = 3.27 
## Functional richness (99.9% probability threshold; Eupjam_GrAn) = 3.98 
## Functional richness (99.9% probability threshold; Euplan_Ma) = 30.32 
## Functional richness (99.9% probability threshold; Eupmin_Ma) = 6.76 
## Functional richness (99.9% probability threshold; Euppec_Ma) = 5.27 
## Functional richness (99.9% probability threshold; Eupruf_Ma) = 5.66 
## Functional richness (99.9% probability threshold; Eupsat_Ma) = 5.77 
## Functional richness (99.9% probability threshold; Euptri_Cois) = 2.72 
## Functional richness (99.9% probability threshold; Euptri_Ma) = 3.18 
## Functional richness (99.9% probability threshold; Eupvio_Cois) = 5.12 
## Functional richness (99.9% probability threshold; Eupvio_Ma) = 13.01
## 
## Functional divergence (Chlcya_Ma) = 0.66 
## Functional divergence (Chlele_Ma) = 0.77 
## Functional divergence (Chlfla_LeAnKa) = 0.39 
## Functional divergence (Chlmus_GrAn) = 0.58 
## Functional divergence (Chlocc_Ma) = 0.37 
## Functional divergence (Chlscl_GrAn) = 0.86 
## Functional divergence (Eupaff_Ma) = 0.6 
## Functional divergence (Eupchl_Ma) = 0.56 
## Functional divergence (Eupful_Ma) = 0.35 
## Functional divergence (Eupgou_Ma) = 0.68 
## Functional divergence (Euphir_Ma) = 0.39 
## Functional divergence (Eupjam_GrAn) = 0.4 
## Functional divergence (Euplan_Ma) = 0.78 
## Functional divergence (Eupmin_Ma) = 0.6 
## Functional divergence (Euppec_Ma) = 0.51 
## Functional divergence (Eupruf_Ma) = 0.6 
## Functional divergence (Eupsat_Ma) = 0.81 
## Functional divergence (Euptri_Cois) = 0.35 
## Functional divergence (Euptri_Ma) = 0.38 
## Functional divergence (Eupvio_Cois) = 0.82 
## Functional divergence (Eupvio_Ma) = 0.83

plot(x = trait.space.pca.eu.group, type = "groups", 
     quant.plot = T, globalContour = T, pnt = T, pnt.cex = 0.1)

Coereba flaveola

coe_traits <- data_traits |>
  filter(Species == "Coereba flaveola")

coe_traits_scl <- as.data.frame(scale(coe_traits[c(3:9)]))

corrplot(round(cor(coe_traits_scl),2), type="upper", order="hclust", 
         tl.col="black", tl.srt=45)

rownames(coe_traits_scl) <- coe_traits$ID

# Calculating the dimensionality of a functional space based on PCA
funspaceDim(coe_traits_scl)

## 
## Using eigendecomposition of correlation matrix.

## 
## Retain the first 1 components

## [1] 1